Exchange and correlation near the nucleus in density functional 

theory 

Zhixin Qian 

Department of Physics, Peking University, Beijing 100871, China 

(Dated: February 6, 2008) 

Abstract 

The near nucleus behavior of the exchange-correlation potential Wxc(r) in Hohenberg-Kohn-Sham 
density functional theory is investigated. It is shown that near the nucleus the linear term of 0(r) of 
the spherically averaged exchange-correlation potential Vxc{r) is nonzero, and that it arises purely 
from the difference between the kinetic energy density at the nucleus of the interacting system and 
the noninteracting Kohn-Sham system. An analytical expression for the linear term is derived. 
Similar results for the exchange Vx{r) and correlation Vc{r) potentials are also obtained separately. 
It is further pointed out that the linear term in fxc(r) arising mainly from fc(r) is rather small, 
and fxc(r) therefore has a nearly quadratic structure near the nucleus. Implications of the results 
for the construction of the Kohn-Sham system are discussed with examples. 

PACS numbers: 71.15.Mb, 31.15.Ew, 71.10.-w 
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Hohenberg-Kohn-Sham density functional theory (HKS-DFT) \^ is now possibly the most 
widely used tool for determining the electronic structure of atoms, molecules, and solids 
In the implementations of HKS-DFT, the so-called Kohn-Sham (KS) noninteracting system 
which has electron density p(r) equivalent to that of the interacting system is assumed. For 
the ground state of an (interacting) A^-electron system in an external potential Vexti^^), the 
one-particle orbitals 0n(i") {n = 1,...N) of the corresponding KS system satisfy the KS 
equation 

[- /2 + Vextir) + vnir) + Vxcir)](l)nir) = e„0„(r), (1) 

with en the eigenvalue. In Eq. ([T]), vnir) is the classical Hartree potential, f_f/(r) = 
J dr'p{r')/\r — r'|, and fxc(r) is the KS exchange-correlation {xc) potential which is the 
derivative of the KS exchange-correlation energy functional E^J^[p] with respect to the den- 

. The KS xc potential takes account of all the cor- 
relations in HKS-DFT. These correlations arise from the Pauli exclusion principle. Coulomb 
repulsion, and correlation-kinetic effects. The former two effects are usually characterized 
with the Fermi-Coulomb hole and the latter arises from the difference between the ki- 
netic energy density of the interacting and the KS systems. E^^^[p] is usually separated into 
two parts, the exchange E^^[p] and correlation Ef^[p] energy functionals. Correspondingly, 
Vrcij^) is split into two components, the KS exchange Vx{r) and correlation fc(r) potentials 

IS, 3- 

In practice, an approximation has to be made for the unknown E^^^[p]. In this respect 
knowledge of the exact properties of fa;c(r) is rather valuable to the construction of the KS 
system. In this paper it is shown that near nucleus, 

Vxc{r) = VxciO) + ^Zim - ts(0)]r/p(0) + ©(r^). (2) 

We use /(r) to denote the spherical average of a function /(r). Here Z is the charge of the 
nucleus, and t{r) and ts(r) are the kinetic energy densities for the interacting system and 
the KS system, respectively. It is further shown that 

VxAr) = + ^Zt,,Mr/p{0) + 0{r^). (3) 

(By the notation Vx,c{''^) we mean Vx{r) and Vc{r) respectively.) Here tx{T) is the first order 
correction to t^fr) in the adiabatic coupling constant perturbation scheme, (see later discus- 
sion and Ref. [5|,) and tc{r) = t^dr) — t^^r). Here txc{r) is defined to be t{r) — ts{r). The 



result in Eqs. ([2]) and ([3]) is valid in general whether the system is an atom, molecule or a 
solid. 

It is well known that,^n the classically forbidden region of finite systems, fa;c(r) decays 
asymptotically as — 1/r [6], a result which has played an important role in the evolution of 
the energy functional construction. In contrast, little is known about Vxc{t^) near the nucleus. 
Inaccuracies in numerical calculations or unreasonable approximations made could even lead 
to divergent results at the nucleus, (see later discussion and Ref. 3]-) It was recently proved 
that fxc(r) is finite at the nucleus. 81] Whether it behaves linearly or quadratically near the 
nucleus is not clear [ol, Q. In this paper it is pointed out that Vxc{t^) approaches the nucleus 
linearly, as shown in Eq. ([2]). 

The asymptotic —1/r structure in the classically forbidden region is totally attributable 



to the Pauli correlation and it has been illustrated to arise from the Fermi hole 



nj. The 



Coulomb correlation only makes a contribution to the next order of 0(l/r^) arising from 



the Coulomb hole 111. It has been shown that the correlation-kinetic effects do make a 

n 

contribution asympotically, but only to the order of 0(l/r^) [11]. Hence in that region the 
correlation-kinetic effects are negligibly small. By contrast Eq. ([2]) indicates that, near the 
nucleus, the linear term of the spherically averaged xc potential v^d''^) arises entirely from 
the correlation-kinetic effects. (This understanding can also be reached via Quantal density 



functional theory (Q-DFT) |l2|. In Ref. 12] it is explicitly shown that the Fermi-Coulomb 
hole makes no contribution to the linear term of Vxc{i^)-) We shall further argue that the 
linear term of fxc(i') is rather small, and fxc(i') is therefore nearly quadratic near the nucleus. 

The understanding reached in this work should be useful in the construction of the KS 
system, especially via approximate xc energy functionals. We shall demonstrate this with 
examples after the main derivations. 

We start with deriving the following cusp relation for the density and the kinetic energy 
density at the nucleus for the interacting system, 

5Zp"(0) + p'"(0) = lz[AZ'piO) + t(0)], (4) 

and the corresponding one for the KS system, 

5Zp"{0) + p"'{0) = ^Z[4ZV(0) + UO) 

+^pmv'H{o) + KM)], (5) 



where the double and triple primes denote the second and third derivatives, respectively, 
with respect to r. To this end, we expand the ground state many-body wavefunction for the 
interacting N electron system as 

^(r, X) = ^(0, X) + a(X)r + b{X.y + c{Xy + ... 
1 

+ [ai^(X)r + hrn{^y]YiUr) + ■■■ 

m=— 1 
2 

+ J2 &2™(X)r2F2„(f) + ..., (6) 

m=-2 

for small r, where f = r/r, and X denotes s, r2S2, • • • , r^SN- By substituting the wavefunc- 
tion into the many-body Schrodinger equation, one obtains the following relations: 

a(X) + Z^(0,X) = 0, 
26i„(X) + Zai^(X) = 0, 

4Z6(X) -Z3^(0,X) + 6c(X) = 0. (7) 

These relations originate from the electron-nucleus cusp of the wavefunction, and were de- 
rived previously in Ref. With them, one can see that the density behaves as 

p(r) = [(1 - Zr)2 + ZVV3]p(0) 

+2(r2 - 5ZrV3)iV j dXi?e[^*(0, X)6(X)] 

m=-l 

+0(r^). (8) 
The kinetic energy density, t(r) = J rfX v X) • V^(i") X), behaves as 

i{r)= m-\zr[2m-z^pm 

-2ZrN j (iXi?e[^*(0,X)6(X)] + 0(r2), (9) 

where 

m = Iz'piO) + N fdXj^ |^l«im(X)p. (10) 

m=—l 

We note that the second term on the right hand side (rhs) of the preceding expression for 
i{0) is absent in the literature j^. Combining Eqs. ([8]) and (fTOj) yields the relation of Eq. 



Next we consider the corresponding KS system. We write 0n(r) as 

0„(r) = ^ [ Anlm + BnlmT + Cni^r^ 
Ira 

+D„z™r3 + ...]l^„(f). (11) 
Substituting Eq. (|TT|) into Eq. ([1]) with fea;t(r) = —Z/r, one can obtain 

-BnOO + ZAnOO = 0, 

'^Bfiiffi + Z A^ijji 0, 

4ZC„oo - {Z' + t)^(O) + t;L(0))A„oo + 6D„oo = 0. 

(12) 

Evidently these relations similarly originate from the electron-nucleus cusp of the KS or- 
bitals. Substituting Eq. (ITT]) into the expression for the density for the KS system, 
P{^) — X]^=i I'/'nl'^)!^) together with the relations in Eq. f[T^ . one has 

p(r) = p(0)[(l - Zrf + \{Z' + v'^iO) + <,(0))r^] 



+^J2R<'^nooCnoo)ir'-^Zr') 

n=l 

1 ^ ' 



n=l m=— 1 

Similarly, the kinetic energy density for the KS system, ts{r) = | V'/'nl'^) ' V0n(r), 

behaves as 

Ur) = isiO) - ^Zr[2t,(0) - Z^p{Q)] 
^ 1 

-2^^^ E 7-^e[A:ooC„oo] + 0(r2), (14) 



Air 

n=l 



where 



n=l m=— 1 

Combining Eqs. (fT3l) and (fTS!) leads to Eq. ([5]). Once again, we note that the second term 
on the rhs of Eq. (fT5l) for ts(0) is absent in the literature 2|. It is exactly this term and 



the corresponding term in t{0) on the rhs of Eq. (fTOj) that make t(0) 7^ ts(0), a fact which 
is critical to the nonzero hnear term of Vxc{r)- In fact, a comparison of Eq. (jl]) and Eq. ([5]) 
leads to t;^(0) + 1;;^(0) = 4Z[t(0) - t,(0)]/3p(0), which equivalently implies Eq. ([2]), since 
v'j,{0) = 0. 

Only the properties of t{0) and ts(0) have been employed in the preceding derivation, 
while accuracy to the linear terms in i{r) in Eq. (Q and is{r) in Eq. (HM enables us to 
obtain an interesting cusp relation, 

t'^M = -^^4c(0), (16) 

with the aid of a comparision of the respective coefficients for the and terms in Eqs. 
dHD and (m. 

To investigate further the near nucleus behavior of the KS exchange Vx{r) and correlation 
tv(r) potentials separately, we employ the adiabatic coupling constant perturbation scheme 
5|. In such a scheme, one has, instead of Eq. ([2]), 

vLir) = vLiO) + lz[t\0) - my/piO) + 0{r'). (17) 



On the other hand, it has been shown that 



1^ 



vi{r) = Xvx{r) + v^{r), (18) 

and v^{r) commences in second order of A. By comparing Eqs. fll7l) and (llSp . one arrives 
at Eq. ([3]). In the meanwhile one has 

= -IzixAo), (19) 

in addition to Eq. ( |T6i) . 

It is well known that tx{r) integrates to zero, i.e., = J drtx{r) = 0. Notice that tx(r) 
is however not necessarily zero except for homogeneous systems. Nevertheless, the fact that 
Tx = somewhat indicates that the linear term of Vxc{r) mainly arises from Vc{r), according 
to Eq. ([3]). In other words, the linear term is rather small, and fa;c(r) is nearly quadratic 
near the nucleus. This argument might be corroborated by the fact that both the following 
sphericalized approximate exchange potentials approach the nucleus quadratically: (a) the 



Pauli correlated approximation Wx{r) of Q-DFT 15|, which is the part of the exchange 
potential arising purely from the Fermi hole, 

Wx{r) = Wx{0) + O{r')- (20) 



(b) the Krieger-Li-Iafrate (KLI) approximation [16] to the optimized exchange potential 



(OEP) (l7|, 



v^''{r)=v^''{0)+O{r'). (21) 
Equation fl20l) follows directly from the fact that the Fermi hole for an electron at the nucleus 

n 

of a sphericalized system is spherically symmetric [12] • Equation (pT]) can be analytically 
derived. (It is worth mentioning that the exact OEP, being the exact exchange-only scheme 

paper.) Therefore, 

in contrast to the common wisdom that Vx{r) dominates fc(r) as is the case in most other 
regions, f c(r) plays a much more significant role in the linear term of fxc(r) near the nucleus. 

A nearly exact result for the density of the helium atom is available, which makes it an 
excellent testing ground for various approximate energy functionals. The single occupied 
KS orbital for the two electrons with opposite spins is simply 0(r) = a/ p(r)/2. The kinetic 
energy density at the nucleus of the KS system can be readily shown, by employing this 
orbital, to be 

tM = Iz'piO). (22) 

Therefore, according to Eq. ([2]), one has f^c(O) = 2Z[2t(0)/p(0) — Z'^]/3. Further employ- 
ment of the cusp relation of Eq. (jll) yields 

<,(0) = [5 Zp"(0) + p"'iO) - 12ZV(0)]/2p(0), (23) 



a result obtained previously in Ref. [10|. Furthermore, in this case, Vx{r) = —VH{r)/2. 
Therefore Vx{r) has no contribution to the linear order, and Vc{r) makes the entire contri- 
bution to the linear term of Vxc{r). Notice that interestingly the second term on the rhs of 
Eq. (|T5l) for ts(0) is absent in Eq. (l22i) in this special case, since both electrons occupy the 
Is orbital. This is in fact also true for three- or four-electron systems like Li and Be atoms 
and their isoelectronic sequences since all electrons occupy single particle states with / = 0. 
The corresponding second term on the rhs of Eq. (fTOj) for t{0) is however nonzero, which 
makes the entire contribution to v'^^{0). 

Table 1 

Comparison of the exact results for Vx^d^) and f^c(O) and those calculated with LDAs for 
the helium atom, (a), (b) and (c) respectively refer to Vosko-Wilk-Nusair, Perdew-Wang, 



and Wigner parametrizations [19] for the correlation energy functionals. Atomic units are 
used. 





VciO) 


<(0) 


<(0) 


Exact -1.69 


-0.062 





3.8 




-0.092 (a) 




0.036 (a) 


LDA -1.51 


-0.082 (b) 


2.02 


0.035 (b) 




-0.055 (c) 




0.0025 (c) 



All the commonly used generalized gradient approximation (GGA) functionals suffer from 



spurious divergence for VxciT") at the nucleus 
GGA in Ref. 







Il0| . This is also true for the later proposed 



18| . The divergence originates from terms in v^d^) containing factor v^p(r) 
which has behavior of 0{l/r) at the nucleus. The local density approximation (LDA), on 
the other hand, yields reasonable result for Vx{0) and Vc{0) as shown in Table 1 for the 
helium atom. But, as shown in Table 1, the LDA for the exchange component largely 
overestimates f^(0), and the three LDAs proposed for the correlation component all hugely 
underestimate f^(0), against the fact revealed in this paper that the correlation component 
plays a dominant role in the linear term of fa;c(r) near the nucleus. 

Finally we take the electrostatic interpretation of the KS xc potential [20^ as another 
example to illustrate the possible implications of our investigation for the construction of 
the KS system. In this interpretation, the concept of static exchange-correlation charge 
density qxd'i') is introduced for the KS xc potential, 

V^Wxc(r) = -47rg^c(r). (24) 

The result shown in Eq. ([2]) indicates that gxc(r) diverges at the nucleus, 

qUr - 0) = --^Z[t(0) - t.(0)]/p(0). (25) 
onr 

This unphysical feature reveals inherent shortcomings of the electrostatic interpretation of 
Wa;c(r). It implies that the part of fxc(r) due to the correlation-kinetic effects can not be 
properly interpreted in terms of a static charge density. Furthermore, according to the result 



in Eq. ([3]), the shortcomings persist in the analogous interpretation 20|] of the exchange ^^^(r) 
and correlation fc(r) potentials, separately, in terms of static exchange ^^(r) and correlation 
gc(r) charge densities. 
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It might be necessary to point out the following fact before we close the paper. As men- 
tioned in the introduction, in addition to the traditional exchange-correlation energy Exc[p\ 
of the conventional many-body theory, which arises purely from the interaction between an 
electron and its Fermi-Colomb hole, there is a contribution of the correlation-kinetic energy 
Txc[p\ = / dr[t{jc) — tsijc)] to the KS exchange-correlation energy: E^f[p\ = E^dp] + T^dp]- 
Correspondingly, Vxdr) = 6Exc[p]/Sp + 6Txc[p]/Sp. The linear term in Vxc{r) in Eq. ([2]) is 
identified as the correlation-kinetic effects, but it remains not clear whether it arises entirely 
from 6Txc[p]/Sp. Nevertheless, the result in this work indicates that, near the nucleus, T2.c[p] 
plays a significant role in the performance of the approximate KS xc energy functionals. 

The author is grateful to Prof. V. Sahni for discussions and help, and Prof. C. J. Umrigar 
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the Chinese National Science Foundation under Grant No. 10474001. 
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